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Abstract 

We demonstrate that two-dimensional (2D) arrays of coupled quantum dots (QDs) with six-fold 
degenerate p orbitals can display bistable states, suitable for application in high-density memory 
device with low power consumption. Due to the inter-dot coupling of px and py orbitals in these 
QD arrays, two dimensional conduction bands can be formed in the x-y plane, while the pz orbitals 
remain localized in the x-y plane such that the inter-dot coupling between them can be neglected. 
We model such systems by taking into account the on-site repulsive interactions between electrons 
in Pz orbitals and the coupling of the localized pz orbitals with the 2D conduction bands formed 
by Px and Py orbitals. The Green's function method within an extended Anderson model is used 
to calculate the tunneling current through the QDs. We find that bistable tunneling current can 
exist for such systems due to the interplay of the on-site Coulomb interactions (U) between the 
Pz orbitals and the delocalized nature of conduction band states derived from the hybridization 
of Px/Py orbitals. This bistable current is not sensitive to the detailed band structure of the two 
dimensional band, but depends critically on the strength of U and the ratio of the left and right 
tunneling rates. The behavior of the electrical bistability can be sustained when the 2D QD array 
reduces to a one-dimensional QD array, indicating the feasibility for high-density packing of these 
bistable nanoscale structures. 
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Intrinsic hysteresis in DC current-voltage characteristics is one of the most intriguing 
problems for resonant tunneling diodes (RTDs)^'^. Such a bistabihty has an important ap- 
plication in memory devices^'^. Whether this phenomenon exists in nanoscale devices such 
as single-electron transistors (SETs) and single molecular transistors (SMTs) has been theo- 
retically investigated in refs [5-7]. Alcxandrov and coworkers^ pointed out that the tunneling 
current through a highly degenerate states of a single QD (molecular) can lead to a switch- 
ing effect only in the case of attractive electron Coulomb interactions, which is mediated 
by electron-phonon interaction. On the basis of Hartree approximation and polaron effect 
Galperin et al proposed that the hysteresis of I-V characteristics can be observed in a single 
molecular junction with effective attractive electron Coulomb interaction.^ Recently, Magna 
and Deretzis showed hysteresis feature of tunneling current in a polaron model beyond the 
Hartree approximation.^ 

Although previous theoretical studies predicted the existence of hysteresis in a QD (or 
molecular) junction, such a phenomenon still lacks conclusive experimental support. 
Moon et al. have experimentally examined the tunneling current through a carbon nanotube 
QD, which exhibits a periodic oscillatory behavior with respect to the applied gate voltage 
arising from the eightfold degenerate state. ^ In addition, Liljeroth et al. have reported a 
periodic oscillatory differential conductance as a result of tunneling current through a single 
spherical PbSe QD with a sixfold degenerate state^. These two experiments did not exhibit 
the bistable tunneling current. Their results indicate that electron-phonon interactions in 
nanotube QDs or PbS QDs are not sufficient to yield the strongly attractive Coulomb inter- 
actions needed for observing the bistabihty. Thus it remains questionable whether a single 
QD junction can display the bistable memory effect. 

Recently, it was demonstrated that semiconductor quantum dot arrays (QDAs) can be 
chemically fabricated to form a superlattice.^°~^^ Via nanoscale manipulation, experimen- 
talists can now control the lattice constant and QD size to tune charges of QDA in the 
Coulomb blockade regime or semiconducting regime. Consequently, QDA is not only 
a good physical system for investigating strongly correlated problem but also a promising 
integrated electronic device. ^'''^^ Although many theoretical efforts have been devoted to the 
charge transport through a single QD,^^'^^ not many studies are on the tunneling current 
through a QDA junction. In this letter we illustrate that a new mechanism exists in a QDA 
junction involving degenerate p-like orbitals which can lead to bistable tunneling current. 
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making it a good candidate for high density storage device. 

Figure. 1 illustrates the system of a QD array embedded in an insulator connected 
with metallic electrodes. The system can be described by the Anderson Hamiltonian, H — 

Hq + Ht + Hd. The Hq = Y,k,a,i3 ^k(i\^a,(3^k,a,p describes the electronic states in the metallic 
leads. Here a\^p (ofc.o-,/?) creates (destroys) an electron of momentum k and spin a with 
energy in the (3 metallic electrode. The Ht term describes the coupling between the 
electrodes and the Pz orbitals of the QD array. 

k,a,P/ k,a,/3/ 

where V^ g,^ describes the coupling between the band states in the electrodes and the localized 
Pz states. Here we assume that the coupling between the electrodes and the Px/Py orbitals of 
the QD array is negligible since the Px/Py orbitals are much more localized along the z axis 
than the Pz orbitals. At last, the term describes electronic states and their interactions 
in the QD array. 

Ha^Y. EpA,A,- + E(^P,A + U{N, - N,))cl,c,,, (2) 
+ Y (^p/4,<tC^V + h.c) + Y UedlJe,adl-adi,-a 

e,p,p',a 

d\ ^ {di^a) creates (destroys) an electron in the pz orbital (with energy Ep^ = Ea) of the 
QD at site £. The second term in Eq. (2) describes the conduction bands of QD array arising 
from the p^ and Py orbitals. A labels the conduction bands (including spin). U denotes the 
on-site Coulomb interaction between two electrons in the p^ and Py orbitals. Note that if we 
ignore the quadrupole and higher-order terms in the expansion of l/ri2, then the Coulomb 
repulsion integrals between two electrons in any of the three degenerate p-like orbitals are 
the same. N\ is the occupation number per unit cell for the A-th conduction band, and 
Nc = J2x is the total occupation number per unit cell for the conduction bands. A mean- 
field theory (which is justified for extended states) has been applied to the 2D conduction 
bands to obtain the second term in the above equation. The third term in Eq. (1) describes 
the hopping coupling between the Pz orbital and the Px/Py orbitals within the tight-binding 
model. The last two terms in Eq. (2) involve Ui — U, and Udc, which denote the on-site 
repulsive Coulomb energy in the Pz orbital, and electron Coulomb interactions between the 
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Pz and Px/Py orbitals. N denotes the number of QDs in the matrix. Here, we focus on 
the Pz orbital rather than the ground state orbital, even though its wave function is more 
localized than that of Pz, since in the range of applied bias considered, the QD ground state 
energy level is deeply below the Fermi levels of both electrodes and the electron tunneling 
through the QD ground state is blockaded. Consequently, carriers in the QD ground states 
only lead a constant-shift to all the p orbitals. It is worth noting that Eq. (2) is similar to 
the so-called the extended Falicov- Kimball model, which has been used extensively to study 
the semiconductor-metal transition in a solid consisting of localized orbitals and delocalized 
orbitals. 

Using Keldysh Green's function technique^^'^^, the tunneling current through the ith QD 
can be expressed by 

Je,. = ^ / -ihie) - /i.(6)] /^'f;;^ ImGl,{e), (3) 

where /l = /(e — i^l) and fu — f{e — /ir) are the Fermi distribution functions for the 
left and right electrodes, respectively. The chemical potential difference between these two 
electrodes is related to the applied bias by /^l — /xr = eVa- r^i(e) and r^fl(e) [F^^^j = 
27rl]k |V£,/3,kp5(e — Ck)] denote the tunneling rates from the pz orbitals to the electrodes. 
Notations e and h denote the electron charge and Plank's constant. In the wide-band limit, 
these tunneling rates are approximately energy-independent. Therefore, the calculation of 
tunneling current is entirely determined by the spectral function A — ImG\^{e), which is 
the imaginary part of the retarded Green's function G\^{e). 
Using the equation of motion for G\p we obtain 



{e-Eo + iV)Gl^{e) 



(4) 



// / 



{e - Eo + ir)Gl^{e) 

= [/ < ni-adi,acl^„ > 

+ E^iy<^p',P+ E 9p",p' <cK^^,Cp>^^,di,acl^>, 

W II I 



(5) 
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(e - e,, A - U{N, - N,))G;,^^{e) (6) 

= Sp',p + Vp>,iGlp+ 9p",p' < KT+^i,i)^^p",a4,a >> 

// 

i,p ,a 

and 

(e-e,,,-C/(7V,-Ar,))G';,^.(e) (7) 

n 

Here, T = (r^,^ + Tg^R)/2 and g^^f = ^^e*(P-P')-^\ In Eqs. (4)-(7), we have intro- 
duced four one-particle Green's functions G[j(e) =< di^^dj^a- ^i,pi^) — < di^T<4),T 
Gp, ,p{e) —< Cp'^acl^a > ^p'j(^) V,CT^]f,(T >■ These four single-particle Green's 
function are coupled with two-particle Green's functions via U and Udc- The equa- 
tion of motion for the two-particle Green's function (defined as < rii-adi^f^dj^^ >, < 
ni,-adi,acl^ >, < cK^y.p'^^,d,^ad]^^ >, < c^^n ydi^acl^ >, < (ni,| + ni^i)dp" ^^cl^ >, and 

< (nj^i + ^i,i)'^p",o-'^j,o- >) coupled to the three-particle Green's functions. In order to 
terminate the heirachy of the equation of motions, we use the Hartree-Fock approximation 
method^^"^^ to decouple terms involving the U^c factor. Meanwhile in the derivation for 

< ni-adi,adj^fj > and < rii-adi^aC^p a >, the treatment for coupling terms between localized 
states and the electrodes (or 2-D conduction band) is employed in the scheme considered 
in our previous method, which is valid for the Coulomb blockade regime.^^'^^ Solving Eqs. 
(4)- (7), we obtain 



€-€p-Ax {e - €p - Ax){€ - €p> - Ax)' 



'^pp'\\^) - z — : — — + - — : — ^-^7: — ^ — 



where = UdA^d^a + ^d,-a) + U{N,- Nx) and 



I — Nd-a 

= ' -;r (9) 

e- Eq- l\c - (rb(e) - iT) 



The retarded Green's function G\^i{e) has the self-energies Ac = U^cNc and Ti,{e) — 
J2p g_g -Ad+iS ^ positive infinitesimal number), which results from the interaction 



between the localized states and conduction band. is the occupation number of Pz or- 
bital in each unit cell. The second term in Eq. (8) describes the scattering amphtude of the 
conduction electron due to interaction with the orbitals. This term would be important 
for studying the charge transport through the Px and Py orbitals in the x-y plane . In this 
study we focus on the longitudinal transport (along z-axis) rather than transverse transport 
(in the x-y plane), thus this term can be ignored. 

To reveal the tunneling current behavior, the occupation number N^^^j determining the 
probability amplitude of resonant channels e — Eq + Ac + (Ti, — iT) and e — Eo + U -|- Ac -|- 
(Ffe — ir) is solved by the equation 



J n i- L + i- R 



As for Nc, we have 



^ _ ^ r ^ r,,j,(e) + r,,j,(.) ^^^^^ ^^^^^^ ^^^^ 

where ^p;^p;,(e) = l/{e- ep^x-U{Nc- Nx)-2UdcNd^i{V l,c + ^ r,c) l'^)- As mentioned above, 
the coupling between the electrodes and the Px/Py orbitals of QD array is negligible, (i.e 
Fc = (ri_c+ri?,c)/2, where ri^c(ri?,c) denotes the tunnehng between the left (right) electrode 
and Px/Py orbitals is small), therefore, IrnQp^p)^{€) fa 7r5(e — €p^x — U{Nc — Nx) — 2UdcNd). 

The range of applied bias considered here would not be enough to overcome the charging 
energy of f/ + Ac, therefore, the second term in Eq. (9) can be ignored and we have 
G^^(e) = (1 — Nd-a) /(e — -Eq ^ Ac — Ff, + iV) in which rb(e) = —iTq (ignoring the small real 
part). The occupation number at zero temperature is calculated by 

^ (i-iv..-,) r. ^^^^ 



Ep+eVa To + r 

de- 



(e-Eo- aeVa - A^)^ + (Fq + Vy 
or 



r 



L 



-7rNd/(l - Nd) (13) 



^_uEf + eVa -Eo- aeVa - UcdNc. 
— cot ( ), 

^ To + r ^' 



in which aeVa term arises from the applied bias crossing QDA and a is a dimensionless 
scahng factor determined by the QDA location, and 

^L,c + ^R,c 



N, = 3^ / deDxie - 2UdcNd - U{Nc - N^)), (14) 



where Dx{e) — Z^p(^(e — ep^\)/N denotes the density of states per unit cell of the A-th 
conduction band. Due to the fact that the Pz energy level is always above the Fermi energy 
of right electrodes (in the range of bias considered), we can ignore the electron injection 
from the right electrode in Eqs. (12) and (14). We first consider the simple case in which 
= (valid for a square lattice) and we approximate the density of states by a square 
pulse function 

D^{e) = Dy{e) = l/W for E^<e<E„ + W, (15) 

where E^, denotes the bottom of the conduction band and W is the band width. Such an 
approximation allows Eq. (11) to have a simple analytic solution of the form 

N^^h- cNa (16) 

with 6 + (1 - a)eVa - Ej]/{-fW + 2,U) and c = 2Udc/{lW + 3[/), where 7 = {Vl,c + 

^R,c)/^L,c- Substituting this into Eq. (13) allows a simple transcendental equation, which 
can be solved numerically. The equation allows a maximum of three roots, out of which 
only two are stable roots. 

We can also solve two coupled transcendental equations as given in Eqs. (13) and (14) nu- 
merically for a more realistic density states, which is derived for a 2D tight-binding model. 
We consider a tight-binding model for and Py orbitals arranged on a rectangular lat- 
tice with lattice constants a and h. Figure 2 illustrates the rectangular lattice. The band 
structure for the Px band is given by 

£^j;(k) — Ep — 2vi cos{kxa) — 2vt cos{kyb), (17) 

where vi denotes the (ppa) interaction and Vt denotes the {pprr) interaction.^^ For the Py 
band, we have 

EyCk) ^ Ep- 2v'i cos{kyb) - 2v[ cos(A;^a). (18) 
The density of states per unit cell form the Px band is given by (if vi > Vt) 

' ^ Jo dri[{2vif - {2vi + 2vt{l - cost]) - if]-^/'^9{i - 2vt{l - cost;)) for < e < Avt 

Dx{e)^l dr][{2vty - {e - 2vt - 2vi{l - cosr))f]-y^e{e - 2vi{l - COST))) for 4vt < e < 4vi , 

^ dr][{2vi)'' - {2vi + 2vt{l + cos 77) - eY]-^/^d{e - 2vt{l + cos 77)) for < e < Avt 

(19) 

where e — e — Ep + 2vi + 2vt and e — Ep + 2vi + 2vt — e. If vi < Vt, then the roles of 
vi and Vt should be exchanged in the above expression. The DOS described by Eq. (19) 



contains the Van Hove singularities. Similar expression {Dy{e)) holds for the Py band with 
the hopping parameters vi and Vt — v\ replaced by v'l — vi and v'^. By varying these hopping 
parameters (for instance, fix lattice constant a and tune b), we can study the behavior of 
bistable tunneling current for systems between the ID and 2D limits. 

We numerically solve the coupled nonlinear Eqs. (13) and (14) for Udc = U — BOmeV, 
Fl — 1 meV{TL,c — Tl/IO), and Tr — ImeViTR^c — Tr/IO). Throughout the paper, we 
shall use T — OX, v[ — 5meV, vi — 20meV, a — 0.5, and £'^ + 1^0 = Ep, where Vq is a, 
reference bias for Va- For simplicity, Vq — 0. The tight-binding parameters are assumed to 
scale according to the rule^^, where R is the separation between two QDs. Thus, we 
have v'l = Viia/h)"^ and Vt = v[(a/by. The occupation number A^^^ as a function of the applied 
bias at zero temperature for the square lattice case (a = b) is shown as solid line in Fig. 3. 
The result remains very similar if we use the constant DOS approximation as described in 
Eqs. (15) and (16) with the same bandwidth, W — 4:{vi + Vt). Although there are Van Hove 
singularities in tight binding DOS, the structure of bias-dependent occupation number does 
not exhibit an anormal feature. This is because Pz orbital is correlated with p^ and Py via 
Nc, which is related to the integral over the DOS. 

Wc sec that the occupation number, has bistable roots. Although, the QDA crystal 
structure reported in references[10,13] has a triangle lattice, the results of Fig. 3 indicate that 
the hysteresis behavior will not depend on the detailed band structure. Once the occupation 
numbers are solved, the tunneling current can be obtained by J = f F^A^d — Jo^d, which is 
valid for zero temperature and when the carrier injection from the right lead can be ignored. 
Consequently, A^^ via the apphed bias directly shows the tunnehng current characteristics. 
The roots for the turn-on and turn-off processes in Fig. 3 are determined by selecting the 
root closest to the root corresponding to the previous value of Va when multiple roots are 
allowed. 

It is crucial to clarify how the system selects the high conductivity state (larger N^) or 
the low conductivity state (smaller No) as the apphed bias is turned on or off. In Fig. 4a, we 
plot the bistable current for various strengths of Udc- Curves 1, 2, and 3 denote, respectively, 
f/=50, 30, and 20 meV. For smaller Coulomb interactions, the bistable current vanishes. 
The critical Coulomb interactions to maintain the bistable current depend on the physical 
parameters such a bandwidth of 2-D conduction band, tunneling rates between dot array 
and electrodes, the broadening Fq, and temperatures. From the application of memory 
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devices, larger U (smaller dot size) is favored because the bistability behavior will be more 
robust against the increase of temperature and broadening. Although Fig. (4a) exhibits the 
bistable current, it does not show any negative differential conductivity (NDC), unlike the 
bistable current in quantum well systems, which is typically associated with NDC.^'^ Fig. 
4b shows the behavior of A^^ for the same set of on-site Coulomb interaction strengths as in 
Fig. 4a. It is useful to understand the behavior of Nc for clarifying the bistable mechanism 
of Nd, which is described below. 

The physical mechanism for this bistablity behavior can be explained as follows. As we 
gradually increase the bias from below the resonance level Ep, the allowed solution to the 
occupation number, A^^^ remains small during the "turn-on" process, as a result of inter-level 
Coulomb blockade (i.e. the level is pushed up by the amount NJJdc-, where the 2-D 
conduction band occupation number is appreciable). Once the applied bias reaches a 
critical value (the inter- level Coulomb blockade is overcome) , A^^ increases quickly to a value 
around 1/3. Now, charges accumulate in the localized Pz orbitals, which leads to an increase 
of the self-energy of the conduction band states by 2UdcNd, and hence reduces Nc to a much 
smaller value (from Ac2 to Acs as illustrated in Fig. 4b). This in turn causes the pz energy 
level to decease during the "turn-off' process while A^^ maintains around 1/3 to keep the 
transcendental equation self-consistent. When the applied bias continues to decrease to a 
critical value, Nc is switched from N^ to Nd (see Fig. 4b), and A^^ quickly goes back to the 
lower value which becomes the only allowed self-consistent solution to the transcendental 
equation. The results of Fig. 4 demonstrate that the bistable tunnehng current arises from 
the on-site repulsive Coulomb interactions Udc- In the polaron model adopted in Refs. 6 and 
7, the bistable mechanism arises from an attractive potential of —-^^d\ ^d£^fjd\ where 
\p and Wq are the electron-phonon interaction strength and the phonon frequency. This 
two particle interaction term should be corrected as {U — ^)dl^^di^^dl _^d£^^^ when the 
intralevel Coulomb interaction is included in a single QD or molecule. Due to the large 
repulsive intralevel Coulomb interaction U, a net attractive electron-electron interactions 
mediated by phonon is difficult to achieve, which may explain why a bistable tunneling 
current through a single QD junction has not been observed. 

To realize nanoscalc memory structures, we need to examine whether this bistable cur- 
rent exists in one dimensional array. Figure. 5 show tunneling current through p^ orbital 
for different ratios of b/a. Other parameters are the same as those for Fig. 3. From two 
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dimensional array to quasi-one dimensional array, the bistability behavior is sustained (al- 
though somewhat weaker in the ID limit). When h/a is greater than 3, (results not shown 
here) the bistable current behaves essentially the same as in the b/a = 3 case. This indicates 
that we already reached the ID limit for b/a ^ 3, and the bistable current still exists. If 
we assume a = 3nm, b = 9nm, and 50 coupled QDs along the chain in the x direction are 
needed to establish a band-like behavior, then the density of the memory device is around 
l/(1350nm^) O.bTB/in^. The operating voltage needed is around lOOmV, which indi- 
cates very low power consumption. Because this is a quantum device, the switching time is 
expected to be comparable to the tunneling rate, which is on the order of 1 THz. 

Finally, wc show in Fig. 6 the tunneling current for various tunneling rate ratios Tl/Th 
for the b/a = 2 case. It is seen that the bistability disappears for F^ = 0.1 meV and 
Tr = 1 meV (shell-tunneling condition). In this case, charges are unable to accumulate 
in the orbitals. Consequently, the effect of Udc is suppressed. On the other hand, for 
Tl — 1 meV and Tr — O.lmeV (shell-fiUing condition), the effect of Udc is enhanced. This 
leads to the wider voltage range of bistable current. From results of Fig. 6, we can control 
the bistable current by adjusting the tunneling rate ratio. 

In conclusion, we have illustrated a novel mechanism for generating the bistable tunneling 
by using a junction involving a 2D periodic array of QDs with six-fold degenerate p-like 
states, which consist of localized orbitals interacting with non-localized p^ and Py orbitals 
via the on-site Coulomb interaction. Due to the interplay of Coulomb blockade effect for 
the locahzed state and the self-energy correction to the 2D conduction bands formed by the 
Px and Py orbitals, a bistable tunneling current with well defined hysteresis behavior can 
be achieved. This bistable current is not sensitive to the details of the band structure, but 
sensitive to the charging energy and the ratio of incoming to outgoing tunneling rate. Such a 
hysteresis behavior arises from a collective effect, not observable in a single QD. It is shown 
that this hysteresis behavior can also exist in one dimensional QDA and very high density 
integrated memory circuits can be realized. 
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Figure Captions 

Fig. 1. Quantum dot array is embedded into an insulated matrix sandwiched between 
two metallic electrodes. 

Fig. 2. Two band (p^ and Py) rectangular lattice with lattice constants a and b. Vi(v[) 
and Vt{v'i) denote, respectively, the electron hopping strength of Px{Py) orbital in a and b. 

Fig. 3. Occupation number as a function of applied bias for Udc — U — bOmeV, Tl — 
1 meV(TL^c — Tl/IO), and F^j = 1 meViVR^c — Tij/lO). Solid line and dashed line denote, 
respectively, tight binding DOS and constant DOS. 

Fig. (4a) Bistable current as a function of applied bias for variation strengths U at zero 
temperature. Fig. (4b) shows N,. for different U . Other parameters are the same those 
of Fig. 3. Note that tunneling current is in units of Jq = eVji/h, which depends on the 
tunneling rate of Tr. 

Fig. 5. Bistable current as a function of applied bias for b/a — 1, 2, and 3. Other 
parameters are the same those for Fig. 3. 

Fig. 6. Bistable current as a function of apphed bias for the various ratios of Fi/Fij and 
b/a = 3. Other parameters are the same as those for Fig. 3. 
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